Compute landcover areas

This notebook computes the landcover areas for each point in 'sample_points' based on a 'landcover' vector map and its defined landcover types. The computed areas are exported to a csv file.

This notebook uses GRASS GIS (7.0.4), and must be run inside of a GRASS environment (start the jupyter notebook server from the GRASS command line).

The code in this notebook was modified from examples in "How to write a Python GRASS GIS 7 addon" developed for a FOSS4G Europe 2015 workshop. The code is contained in the python-grass-addon repository on GitHub.com.

Required packages

Variable definitions

landcover – vector map representing land cover
sample_points – vector map containing observer points
radius – radius in map units around each sample point to calculate land cover areas within
use_viewshed – consider only 'visible' land cover, based on a viewshed from each sample point
landcover_area_table_filename — path to export land cover area table as a csv file
n_types — number of land cover types


In [ ]:
landcover = 'landcover'

In [ ]:
landcover_type_field = 'type'

In [ ]:
n_types = 15

In [ ]:
sample_points = 'sample_points_field'

In [ ]:
radius = 500

In [ ]:
use_viewshed = True

get viewshed suffix (for filenames)


In [ ]:
def getViewshedSuffix():
    if use_viewshed:
        viewshed_suffix = '_viewshed'
    else:
        viewshed_suffix = ''
    return viewshed_suffix

In [ ]:
landcover_area_table_filename = "/home/ubuntu/landcover_areas_{0}m{1}.csv".format(radius, getViewshedSuffix())

Import statements


In [ ]:
import numpy as np
import pandas
import pyprind

GRASS import statements


In [ ]:
import grass.script as gscript

In [ ]:
from grass.pygrass.vector.geometry import Point
from grass.pygrass.vector import Vector

from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.table import DBlinks

Function declarations

connect to attribute table


In [ ]:
def connectToAttributeTable(map):
    vector = VectorTopo(map)
    vector.open(mode='r')
    dblinks = DBlinks(vector.c_mapinfo)
    link = dblinks[0]
    return link.table()

extract point from vector map


In [ ]:
def extractPoint(input, ID, output):

    where = 'ID = {0}'.format(ID)
    type = 'point'

    gscript.read_command('v.extract',
                         input=input,
                         where=where,
                         output=output,
                         type=type, 
                         overwrite=True)

create buffer around point


In [ ]:
def bufferPoint(input, output, radius):
    gscript.read_command('v.buffer',
                         input=input,
                         output=output,
                         type='point',
                         distance=radius,
                        overwrite=True)

convert viewshed from raster to vector map


In [ ]:
def vectorizeViewshed(input, ID, output):
    type = 'area'
    column = 'visible'

    gscript.read_command('r.to.vect', 
                         input=input, 
                         output=output,
                         type=type,
                         column=column,
                         overwrite=True)

overlay a vector map on an underlying vector map using 'and' selection operator


In [ ]:
def overlay(overlay, underlay, output):
    operator='and'
    gscript.read_command('v.overlay',
                         ainput=overlay,
                         binput=underlay,
                         operator=operator,
                         output=output,
                         overwrite=True)

add area column to vector map


In [ ]:
def calculateAreas(map):
    
    #add new area column
    gscript.read_command('v.db.addcolumn',
                         map=map,
                         columns="area_square_meters DOUBLE PRECISION")
    
    #compute area and insert into area column
    gscript.read_command('v.to.db',
                         map=map,
                         type='centroid',
                         option='area',
                         columns='area_square_meters',
                         unit='meters')

create table showing total area by landcover type


In [ ]:
def getLandcoverAreaByType(map):
    #get area data
    table = connectToAttributeTable(map=map)
    table.filters.select()
    columns = table.columns.names()
    cursor = table.execute()
    result = np.array(cursor.fetchall())
    cursor.close()
    data = pandas.DataFrame(result, columns=columns).set_index('cat')

    #sum areas by landcover type
    data['area_square_meters'] = pandas.to_numeric(data['area_square_meters'])
    areas = data[['b_' + landcover_type_field, 'area_square_meters']].groupby(by='b_' + landcover_type_field).sum()
    return areas

Utility functions

export vector map to a shapefile


In [ ]:
def exportVectorToShapefile(map, output):
    gscript.read_command('v.out.ogr',
                         input=map,
                         format='ESRI_Shapefile',
                         output=output,
                         flags='e',
                         overwrite=True)

get info about a vector map


In [ ]:
def getVectorMapInfo(map):
    return gscript.read_command('v.info', map=map)

In [ ]:
#view(['1_viewshed'])
# TODO determine how to display vector map

Compute landcover areas

connect to 'sample_points' attribute table


In [ ]:
point_table = connectToAttributeTable(map=sample_points)
point_table.filters.select()
columns = point_table.columns.names()
cursor = point_table.execute()
result = np.array(cursor.fetchall())
cursor.close()
point_data = pandas.DataFrame(result, columns=columns).set_index('cat')

loop through sample points


In [ ]:
with Vector(sample_points, mode='r') as points:
    
    #setup progress bar
    progress_bar = pyprind.ProgBar(points.n_lines, bar_char='█', title='Landcover analysis progress', monitor=True, stream=1, width=50)
    
    #iterate through points
    for point in points:
        
        #get point ID (SiteID)
        ID = point_data['ID'][point.cat-1]
        
        #update progress bar
        progress_bar.update(item_id=ID)
        
        #buffer current point
        extractPoint(input='sample_points_field', ID=ID, output='tmp_buffer_point')
        bufferPoint(input='tmp_buffer_point', output='tmp_point_buffer', radius=radius)
        
        #set buffer as overlay
        overlay_input = 'tmp_point_buffer'
        #consider only visible landcover if 'use_viewshed' = True
        if use_viewshed:
            viewshed = 'vect_{0}_viewshed'.format(ID)
            visible_viewshed = 'vect_{0}_viewshed_{1}m'.format(ID, radius)
            #vectorize viewshed
            vectorizeViewshed(input='{0}_viewshed'.format(ID), ID=ID, output=viewshed)
            #overlay buffer on viewshed
            overlay(overlay=overlay_input,
                    underlay=viewshed,
                    output=visible_viewshed)
            #set overlay to the visible viewshed
            overlay_input = visible_viewshed
        overlay_output = 'vect_{0}_landcover_{1}m{2}'.format(ID, radius, getViewshedSuffix())
        #overlay landcover
        overlay(overlay=overlay_input,
                         underlay=landcover,
                         output=overlay_output)
        
        #calculate landcover area
        calculateAreas(map=overlay_output)

Export landcover area table

create table with landcover areas by type for each point


In [ ]:
#create table
index_start = 41
''' set the first index number,
    allowing easier insertion into a database table that already contains 
    area calculations with other parameters
    (i.e radius and use_viewshed)'''
columns = ['ID', 'SiteID', 'IncludedArea']
columns = columns + [ str(n) for n in range(1,n_types+1) ]
area_table = pandas.DataFrame(columns=columns)

#set naming variables
included_area = '{0}m{1}'.format(radius, getViewshedSuffix())  #0=radius, 1=viewshed_suffix
map_pattern = 'vect_{0}_landcover_{1}m{2}'  #0=ID, 1=radius, 2=viewshed_suffix

#iterate through points
for index, point in point_data.iterrows():
    ID = point['ID']
    map = map_pattern.format(ID, radius, getViewshedSuffix())
    #initiate row
    row = {'ID':"{0:.3g}".format(int(index) + index_start),
           'SiteID': str(ID),
           'IncludedArea': included_area}
    #get landcover areas
    areas = getLandcoverAreaByType(map)
    
    #iterate through area types
    for index, area in areas.iterrows():
        #add area to row
        row["{0:.3g}".format(int(index))] = area['area_square_meters']
    #append row to table
    area_table = area_table.append(row, ignore_index=True)

area_table.set_index('ID', inplace=True)

#export table to file
area_table.to_csv(landcover_area_table_filename, header=False)